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Abstract 

Experimentalists now measure intense rotations of Lagrangian particles in turbulent flows 
by tracking their trajectories and Lagrangian-average velocity gradients at high Reynolds 
numbers. This paper formulates the dynamics of an orthonormal frame attached to each 
Lagrangian fluid particle undergoing three-axis rotations, by using quaternions in combina- 
tion with Ertel's theorem for frozen-in vorticity. The method is applicable to a wide range 
of Lagrangian flows including the three-dimensional Eulcr equations and its variants such as 
ideal MHD. The applicability of the quaterionic frame description to Lagrangian averaged 
velocity gradient dynamics is also demonstrated. 
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1 Introduction 

1.1 General background 

On 16th October, 1843, Hamilton wrote the expression for quaternion multiplication, 

I 2 = f = K 2 = UK = -Id , 

as an algebraic rule for composing "turns" (directed arc lengths of great circles) in orienting his 
telescope. This feature - that multiplication of quaternions represents composition of rotations 
- has made them the technical foundation of modern inertial guidance systems in the aerospace 
industry where tracking the paths and the orientation of satellites and aircraft is of paramount 
importance [1]. The graphics community also uses them to control the orientation of tumbling 
objects in computer animations because they avoid the difficulties incurred at the north and 
south poles when Euler angles are used [2]. 

Given the utility of quaternions in explaining the dynamics of rotating objects in flight, one 
might ask whether they would also be useful in tracking the orientation and angular velocity of 
Lagrangian particles in fluid dynamical situations. Recent experiments in turbulent flows have 
developed to the stage where the trajectories of tracer particles can be detected at high Reynolds 
numbers [3, 4, 5, 6, 7, 8, 9, 10, 11, 12] ; see Figure 1 in [3]. Numerical differentiation of these 
trajectories gives information about the Lagrangian velocity and acceleration of the particles. 
In particular the curvature of the particle paths can be used to extract statistical information 
about velocity gradients from a single trajectory [12]. 

The usual practice in graphics problems is to consider the Frenet-frame of a trajectory which 
consists of the unit tangent vector, a normal and a bi-normal [2, 12]. In navigational language, 
this represents the corkscrew-like pitch, yaw and roll of the motion. In turn, the tangent vector 
and normals are related to the curvature and torsion. While the Frenet-frame describes the path, 
it ignores the dynamics that generates the motion. Here we will discuss another ortho-normal 
frame associated with the motion of each Lagrangian fluid particle, designated the quaternion- 
frame. This may be envisioned as moving with the Lagrangian particles, but their evolution 
derives from the Eulerian equations of motion. 

The first main contribution of this paper lies in the explicit calculation of the evolution 
equations of this quaternion frame. Secondly, this formulation is shown to apply to a wide range 
of Lagrangian problems, as the next sub-section shows. Thirdly, the evolution of the quaternion 
frame is related to an associated Frenet frame. Finally, it is shown how the pressure Hessian, 
which plays a key role in driving the quaternion frame for Euler fluid flow, can be modelled 
using constitutive relations. 



1.2 An appropriate class of Lagrangian flows 

Suppose id is a contravariant vector quantity attached to a particle following a flow along 

characteristic paths dx/alt = u(x, t) of a velocity u. Let us consider the abstract Lagrangian 
flow equation 

Dw . , D d _ .„ „. 

= a(x,t), — = — +u-V. (1.1) 

Dt y ' Dt dt y ' 

where the material derivative has its standard definition. Examples of systems that (1.1) might 

represent are : 
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1. The vector w in (1.1) could be the velocity of a tracer particle in a fluid transported by a 
background velocity field u with a as the particle's acceleration. 

2. w could be the vorticity u> = curlu of the incompressible Euler fluid equations, in which 
case a = ld ■ Vit and divu = 0. The case of the Euler equations with rotation ft would 
make w = u = p~ 1 (uj + 2 ft). 

3. For the barotropic compressible Euler fluid equations (where the pressure p = p(p) is only 
density dependent) then w = u p = p~ x uj, in which case a = u p ■ Vu and divu =/ 0. 

4. The vector w could represent a small vectorial line element S£ transported by a background 
flow u, in which case a = S£ ■ Via. For example, using Moffatt's analogy between the 
(magnetic) S-field in ideal MHD and fluid vorticity [13], taking S£ = B in the equations 
for incompressible ideal MHD, then a = B ■ Vu and div_B = 0. In a slightly generalized 
form it could also represent the Elsasser variables = u±B in which case a ± = w^ Vu 
with two material derivatives [13]. In each of the cases (2-4) the vectors w and u satisfy 
the standard Eulerian form 

Dw „ 

= w-Vu. (1.2) 

Dt v ; 

Consequently, it follows from Ertel's Theorem [14, 15, 16] that 

for any differentiable function /i(cc, t). Choosing fj, = u(x, t) as in [15] and identifying the 
flow acceleration as Du/Dt = Q(x, t) yields the second derivative relation 

D 2 w (Du\ „ , . 

w-V ( — ) =:w-VQ. (1.4) 



Dt 2 V Dt , 

In each of the cases (2-4) above the acceleration vector Q in (1.4) is readily identifiable. Thus, 
in those cases we have 

Dw Da Dx 

-j^j- = a(x, t) and = w ■ \/Q =: b(x, t) along = u(x, t) . (1-5) 

These are the kinematic rates of change of the vectors Dw/Dt = a and Da/Dt = b following the 
characteristics of the velocity vector u along the path x(t) determined from Dx/Dt = u(x, t). 

The plan of the paper is as follows: §3 shows that the quartet of 3- vectors (u, w, a, b) 
appearing in (1.5) determines the quaternion-frame and its Lagrangian dynamics: indeed, a 
knowledge of b is essential to determining the dynamical process. Modulo a rotation around 
w, the quaternion- frame turns out to be the Frenet-frame attached to lines of constant w. In 
case (2) where w = u>, lines of constant w are vortex lines and the particles are fluid parcels. 
As described in [15] - see also [16] for a history - in this case Ertel's Theorem for Euler's fluid 
equations [14] ensures that b = — Pu where P is the Hessian matrix of spatial derivatives of the 
pressure. In some practical situations, however, the vector b may not be known, or may not 
exist for every system for every triad (u, w, a). For example no explicit relation for b is known 
for the Euler equations in velocity form for which (u, w, a) = (u, u, — Vp). 

Section §4 elaborates three examples: the Euler equations with rotation; the barotropic 
compressible Euler equations; and ideal MHD in Elsasser variables. Section §5 demonstrates the 
applicability of the quaterionic frame description to turbulence models of Lagrangian averaged 
velocity gradient dynamics, in which approximate constitutive relations for the pressure Hessian 
and auxiliary equations for the constitutive parameters are introduced. 



16/02/07 



http://arxiv.org/abs/nlin.CD/0607020 



4 



2 Quaternions and rigid body dynamics 

Rotations in rigid body mechanics has given rise to a rich and long-standing literature in which 
Whittaker's book is a classic example [17]. This gives explicit formulae relating the Euler angles 
and what are called the Cayley-Klein parameters of a rotation. In fact the use of quaternions in 
this area is not only much more efficient but avoids the immensely complicated inter-relations 
that are unavoidable when Euler angle formulae are involved [18, 19]. 

In terms of any scalar p and any 3-vector q, the quaternion q = [p, q] is defined as (Gothic 
fonts denote quaternions) 

3 

q = \p, q] =p\ -J^fto-i, (2.1) 

i=l 

where {o~x, 02, 03} are the three Pauli spin-matrices defined by 

-(!;)• -UO" m 

and I is the 2x2 unit matrix. The relations between the Pauli matrices <jj<jj = — <%l — tijk&k 
then give a non-commutative multiplication rule 

qi © q2 = [P1P2 - <?i • <?2, Pi<?2 + P29i + Qi x qr 2 ] • (2.3) 

It can easily be demonstrated that quaternions are associative. 

Let p = \p, q] be a unit quaternion with inverse p* = [p, — qr]: this requires p © p* = 
[p 2 + q 2 , 0] = [1,0] for which we need p 2 + q 2 = 1. For a pure quaternion r = [0, r] there exists 
a transformation from r = [0, r] — ► 9\ = [0, i£] 

9r = p©t©p*. (2.4) 

This associative product can explicitly be written as 

JR = p © r © p* = [0, (p 2 - q 2 )r + 2p(q x r) + 2g(r • q)} . (2.5) 
Choosing p = ± cos |0 and q = ±n sin where n is the unit normal to r, we find that 

$H = p © t © p* = [0, r cos 9 + (n x r) sin 0] = O(0 , h)r , (2.6) 

where 

p = ±[cosi0, n sin 10] . (2.7) 

Equation (2.6) is the famous Euler- Rodrigues formula for the rotation O(0 , n) by an angle 
of the 3-vector r about its normal h ; the quantities , n are called the Euler parameters. The 
elements of the unit quaternion p are the Cayley-Klein parameters 1 which are related to the 
Euler angles [17]. 

Lemma 1 The unit quaternions form, a representation of the Lie group SU(2). 



The Cayley-Klein parameters of the quaternion q = [a, x] of §3 are given by 
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Proof : From (2.1), the matrix representation of a unit quaternion is 

( p-iqz -iqi - q 2 \ , n ^ 

J = p I — iq ■ cr = (2.8) 
y -zgi + q 2 p + iqz J 

where J £ SU(2) ; that is, J is a unitary 2x2 matrix with unit determinant. Hence, we may 
rewrite the map (2.4) for quaternionic conjugation equivalently in terms of the Hermitian Pauli 
spin matrices as 

R a = Jr crjt (2.9) 

This is the standard representation of 50(3) rotations as a double covering (±J) by 577(2) 
matrices, which is now seen to be equivalent to quaternionic multiplication : for more discussion 
of this theorem see a modern treatise on mechanics such as [20]. The (±) in the Cayley-Klein 
parameters reflects the 2:1 covering of the map 577(2) i— ► 50(3). ■ 

To investigate the map (2.4) when p is time-dependent, the Euler-Rodrigues formula in (2.6) 
can be written as 

9\(t)=p®x®p* x = p* ©9t(i) ©p . (2.10) 

Thus yi has a time derivative given by 

fk(t) = p ® (p* ©9\©p) ©p* + p ® (p* ®<R®p) ®p* 

= p©p*©9\ + 9l©p©p* 

= (p ©p*) ©9i + ^© (p ©p*)* 

= (^©p*)©^-((^©P*)©^)*, (2.H) 

having used the fact on the last line that because 91 is a pure quaternion, 9\* = — 9t. Because 
p = [p, q] is of unit length, and thus pp + qq = 0, this means that p©p* is also a pure quaternion 

£®p* = [0, in (t)]. (2.12) 

The 3-vector entry in (2.12) defines the angular frequency Qo(t) as CIq = 2(—pq + qp — q x q) 
thereby giving the well-known formula for the rotation of a rigid body 

R = £l xR. (2.13) 

For a Lagrangian particle, the equivalent of f2o is the Darboux vector T> a in Theorem 1 of §3. 

3 Lagrangian evolution equations and the quaternion picture 
3.1 An ortho-normal frame and particle trajectories 

Having set the scene in §2 by describing some of the essential properties of quaternions, it is now 
time to apply them to the Lagrangian relation (1.1) between the two vectors w and a which we 
shall repeat here 

Dw 

— = a. (3.1) 

Through the multiplication rule in (2.3) quaternions appear in the decomposition of the 3- 
vector a into parts parallel and perpendicular to another vector, which we choose to be w. This 
decomposition is expressed as 

a = a a w + Xa xw = [a a , x a ] ® [°> w ] > ( 3 - 2 ) 
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where the scalar a a and 3-vector x a are defined as 



a n = w (w ■ a) 



X a = w {wxa) 



(3.3) 



Equation (3.2) thus shows that the quaternionic product is summoned in a natural manner. It 
is now easily seen that a a is the growth rate of the scalar magnitude (w = \w\) which obeys 



Dw 
~Dt 



while x a i the swing rate of the unit tangent vector w = ww , satisfies 

Dw 



Dt 



Xa XW - 



(3.4) 



(3.5) 




W X Xa 



(X2,t 2 ) 



tracer particle trajectory 




Figure 1: The dotted line represents the tracer particle (•) path moving from (xi,ti) to (x2,i2)- 
The solid curves represent lines of constant w to which to is a unit tangent vector. The orientation of 
the quaternion-frame (w, x a , w x x a ) is shown at the two space-time points; note that this is not the 
Frenet-frame corresponding to the particle path but to lines of constant w. 



Now define the two quaternions 

q a = [a a , X a ] 



to = [0, w] , 



(3.6) 



where to is a pure quaternion. Then (3.1) can automatically be re- written equivalently in the 
quaternion form 



Dm 

~Dt 



q a © W 



Moreover, if a is differentiable in the Lagrangian sense (see (1.5)) 

Da 

m =b > 

then, exactly as for q a , a quaternion q can be defined which is based on the variables 
a b = w~ l {w ■ b) , x 6 = tiT 1 ^ x fa) , 

where 

qb = K, Xb\ ■ 

It is now clear that there exists a similar decomposition for b as that for a as in (3.2) 

D 2 K> 



Dt 2 



[0, b] = [0, otbW + Xb x w ] = qb © tt> • 



(3.7) 

(3.8) 

(3.9) 
(3.10) 

(3.11) 
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Using the associativity property, compatibility of (3.11) and (3.7) implies that (w = \w\ /= 0) 

Dc\ a 



Dt 

which establishes a Riccati relation between q a and q b 

Dq a 



+ q ®q -q& ®» = 0, (3.12) 



m +c\a®q a = q b , (3.13) 
whose components yield 

^[«a, Xa] + [«a ~ xl, ^aXa] = K , Xb] , (3.14) 

where Xa = \x a \- Fro 111 (3-13), or equivalently (3.14), there follows the first main result of the 
paper: 

Theorem 1 The ortho-normal quaternion-frame (w, Xai wxx a ) G SO (3) has Lagrangian time 
derivatives expressed as 

Dw 

-- T> a xw, (3.15) 



Dt 

D{W X Xa) 

Dt 

DXa 



-Da^(wX Xa ), (3-16) 
Va^Xa, (3-17) 



Dt 

where the Darboux angular velocity vector T> a is defined as 

?>a = Xa + —W , C b = W • (X a X Xb) ■ (3-18) 

Xa 

Remark : The Darboux vector T> a sits in a two-dimensional plane and is driven by the vector 
b which sits in Q, in (3.18). The analogy with rigid body rotation expressed in (2.13) is clear. 
Proof : To find an expression for the Lagrangian time derivatives of the components of the 
frame (w, Xai x Xa) requires the derivative of x a - To find this it is necessary to use the fact 
that the 3-vector b can be expressed in this ortho-normal frame as the linear combination 

vj-^-b = a b w + c b x a + d b (w x x a ) ■ (3-19) 

where c b is defined in (3.18) and d b = — (x a • Xb)- The 3-vector product x b = ^ _1 (i& X b) yields 

X b = c b (w x Xa) -dbXa- (3-20) 

When split into components, equation (3.14) becomes 



Dctg 
Dt 

and 

DXa . 

Dt 

Prom the latter it is easily seen that 

Dt 



X 2 a -a 2 a + a b (3.21) 
-ZUaXa + Xb ■ (3.22) 

-2a a Xa - d b (3.23) 
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from which it follows 

^ = C b XaHw X Xa) , D{W p t =XaW~ C^Xa , (3-24) 

which gives equations (3.15)-(3.18). ■ 

Theorem 1 is the main result of the paper and is the equivalent for a Lagrangian particle 
undergoing fluid motion of the well-known formula (2.13) for a rigid body undergoing rotation 
about its center of mass. 

3.2 The evolution of the 6-field in equation (3.8) 

The Lagrangian rate of change of acceleration Da/Dt = b is important for tracking passive 
tracer particles. However, the vector b cannot be calculated directly from Ertel's Theorem in 
(1.5). As shown below, the Lagrangian evolution of q b appearing in the quaternionic Riccati 
relation (3.13) may be described without approximation in terms of three arbitrary scalars. 

Theorem 2 The Lagrangian time derivative of the quaternion q b in the Riccati relation (3.13) 
can be expressed as 

^ = q« © q fe + %a,b , (3.25) 

< Pa,fe = Aiq 6 + A 2 q a + A 3 Id, (3.26) 

where X\(x, t) , \2(x, t) , Xs(x, t) are arbitrary scalar functions and Id = [1, 0] is the identity 
for the quaternions. 

Remark: Without further constraints \\{x, t), \2{x, t) and \%{x, t) would be arbitrary. 

Proof: To establish (3.25), we differentiate the orthogonality relation % b • it) = and use the 
Lagrangian derivative of w 



Dt 



- = Xa x Xb + «o , where s = fj, Xa + ^Xb- (3-27) 



so lies in the plane perpendicular to w in which x a an d Xb a ^ so ne an d A* = ^( x : t) an d A = \(x, t) 
are arbitrary scalars. Explicitly differentiating Xb = u] (w x b) gives 

w^w (X a -b) + s = -a a x b - a b x a + w^w (x a ■ b) + uT 1 x , (3.28) 

which can easily be manipulated into 

w x |^ - a b a - a a fo| = w s . (3.29) 

This means that 

Db . 

-p- = a b a + a a b + so x w + ew , (3.30) 

where e = e(x, t) is a third unknown scalar in addition to \i and A in (3.27). Thus the Lagrangian 
derivative of a b = w~ x {w ■ b) is 

Da b 

-j^- = a a a b + Xa ■ Xb + e - (3-31) 
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Lagrangian differential relations have now been found for Xb an< ^ a b, but at the price of intro- 
ducing the triplet of unknown coefficients fj,, A, and e which are re-defined as 

A = a a + Ai , n = a b + A 2 , e = -2x a • Xb + ^a a + Aia 6 + A 3 . (3.32) 

The new triplet has been subsumed into the tetrad defined in (3.26). Then (3.27) and (3.31) 
can again be written in the quaternion form (3.25). ■ 

In §5 we shall discuss a approach for determining by introducing an approximate constitu- 
tive relation for the pressure Hessian P and auxiliary equations for the constitutive parameters. 



3.3 Frame dynamics and the Frenet equations 




Figure 2: The ortho-normal frame (w, w x x a , % a ) as the Frenet- frame to lines of constant w. 

Modulo a rotation around the unit tangent vector w, with x a as the unit bi-normal b and w x x a 
as the unit principal normal n, the matrix F can be formed 



F = (w T , (to x Xa ) T , Xa) 
and (3.15)-(3.17) can be re-written as 

/ 



(3.33) 



151 



AF. 



A 







Xa 

X -CbXa 
\ C^a 1 



(3.34) 



For a space curve parameterized by arc-length s, then the Frenet equations relating dF/ds to 
the curvature k and the torsion r of the line of constant w are 



ds 



where 



N 






K 





— K 





T 





— T 






(3.35) 



It is now possible to relate the t and s derivatives of F given in (3.34) and (3.35). At any time t 
the integral curves of the vorticity vector field define a space-curve through each point x. The 
arc- length derivative d/ds is defined by 



(3.36) 
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The evolution of the curvature k and torsion r of a vortex line may be obtained from Ertel's 
theorem in (1.3), expressed as the commutation of operators 



d D 

ds' Dt 



d 

a a — . (3.37) 

ds 



Applying this to F and using the relations (3.34) and (3.35) establishes the following Theorem 

Theorem 3 The matrices N and A satisfy the Lax equation that relates the evolution of the 
curvature k and the torsion r to a a , Xa and Q, defined in equations (3.2) and (3.18) 

^-a a N = ^ + [A,N]. (3.38) 
Thus, if x a = the curvature k is stationary. 

4 Three further examples 

The quaternionic formulation can be applied to other situations, such as the stretching of fluid 
line-elements, incompressible and compressible motion of Euler fluids and ideal MHD [21]. 

(i) The Euler fluid equations for incompressible flow in a rotating frame: The velocity 
form of Euler's equations for an incompressible fluid flow in a frame rotating at frequency ft is 

Du 

— = (u x 2ft) -Vp , with divu = 0. (4.1) 

Coriolis 

Taking the curl yields (u = curl it) 

— = u ■ Vu , with u = p~ x {u + 2fi) . (4.2) 

The triad of vectors (u, w, a) in this case represents (it, to, u ■ Vu) with u = curlu and 
div u = 0. Then Ertel's theorem becomes 

£< a -™- a -*(^)- (43> 

Upon taking p, = u in Ertel's theorem as in [15] and using the motion equation gives a relation 
in the moving frame 

-^■(cj-Vu) = Q ■ V(« x 2ft - Vp) 

= -Pu> + u ■ V(u x 20) , (4.4) 

where P is the Hessian matrix of the pressure defined by 

d 2 p 



dxidxj 



(4.5) 



Thus (4.3) and (4.4) identify a and b as a = Cb • Vu and b = — P<D + Cb ■ V [u x 2fi) ; the particles 
are now fluid packets and not passive tracer particles. The divergence-free constraint divu = 
implies that 

-Ap = u id u jt i - div (u x 2fi) = tr S 2 - \u? - div (u x 2ft) . (4.6) 
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Equation (4.6) places an implicit condition upon the relation between the strain matrix S and the 
pressure Hessian P in addition to the Riccati equation (3.12). This situation has been discussed 
at greater length in [16] in the absence of rotation. 

(ii) Euler's equations for a barotropic compressible fluid: The pressure of a barotropic 
compressible fluid is a function of its mass density p, so it satisfies V p x Vp = 0. The velocity 
form of Euler's equations for incompressible fluid motion in a frame rotating at frequency ft is 

^ = --Vp(p)=:-Vh(p), (4.7) 

^ + pdivu = 0. (4.8) 

Taking the curl yields 

- = uj p -Vu, with Up = <jj/p and cj = curlti. (4-9) 

Then Ertel's theorem takes the same form as above, and the second Lagrangian time derivative 
yields the Ohkitani relation for a barotropic compressible fluid, 

= ^ P 'V«) = - Wp 'V (V%)) = b , (4.10) 

in terms of the Hessian of its specific enthaply, h{p). This has the same form as for incompressible 
fluids, except the acceleration term b = — u p ■ V (V/i(/o)) has its own dynamical equation. Thus, 
the methods in [16, 22] also apply for barotropic fluids. For isentropic compressible fluids, the 
situation is more complicated. 

(iii) The equations of incompressible ideal MHD: These are 

^ = B-VB-Vp, (4.11) 

^ = B-Vu, (4.12) 

together with divu = and div£? = 0. The pressure p in (4.11) is p = pj + \B 2 where pj is 
the fluid pressure. Elsasser variables are defined by the combination 

v ± = u±B. (4.13) 

The existence of two velocities v ± means that there are two material derivatives 

= ^- + v ± -V. (4.14) 

Dt at v ' 

In terms of these, (4.11) and (4.12) can be rewritten as 

Dt ~ 

with the magnetic field B satisfying (divv^ = 0) 

D ± B 



Vp, (4.15) 



Dt 



B-Vv ± . (4.16) 
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Thus we have a pair of triads (v^, B, a^) with a ± = B ■ \7v^, based on Moffatt's identification 
of the S-field as the important stretching element [13]. From [16, 22] we also have 

^ = -PB, (4.17) 

where b^ 1 = —PB. With two quartets (v^, B, , b), the results of Section 2 follow, with two 
Lagrangian derivatives and two Riccati equations 

^ + q a ± ©q? = q fe . (4.18) 

In consequence, MHD-quaternion-frame dynamics needs to be interpreted in terms of two sets of 
ortho-normal frames (^B, B x x^ acted on by their opposite Lagrangian time derivatives. 

D T B 

= V^xB, (4.19) 

Dt 

1)T (Bx X ± ) = V^xiBxx*), (4.20) 



Dt 



Dt 

where the pair of Elsasser Darboux vectors T)^ are defined as 



V = -D^xx^, (4.21) 



V T = X T ~ %B , 4 = B-[ X ± x (x pB + « ± X T )] • (4-22) 



5 Approximate constitutive relations for the pressure Hessian 

In this section, we discuss an approach for determining q& by introducing an approximate con- 
stitutive relation for the pressure Hessian P in the Euler equations and auxiliary equations for 
these constitutive parameters. Recall Euler's familiar equations for an incompressible fluid flow 
with velocity u, written as 

Du 

-j^j- = — Vp , with divu = 0. (5-1) 

Taking the gradient yields the matrix Riccati equation 

DM 



Dt 



+ P + |VH = 0, (5.2) 



where the velocity gradient tensor M = Vu has Cartesian components Mij = duj/dx 1 = ujj 
and the (symmetric) pressure Hessian P = VVp has Cartesian components Pjj = d 2 p/ dx l dxK 

Because of the incompressibility condition div u = 0, the trace of the velocity gradient tensor 
tr M vanishes, thereby requiring tr P = — tr(M 2 ), which is a Poisson equation for the pressure. 
For laminar flow in a bounded domain, the Poisson equation determines both the pressure in 
the exact Euler equations and its Hessian appearing in the velocity gradient equations (5.2). 

However, in turbulent flows, modern diagnostics for both numerical simulations and fluid 
measurements make extensive use of average values of the velocity gradients moving in a coarse- 
grained volume element following the mean flow [23]. The process of averaging following a 
fluid parcel is called Lagrangian averaging. By its definition, Lagrangian averaging commutes 
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with the material derivative, but with not the spatial gradient. In contrast, Eulerian averaging 
does the opposite. With its two Eulerian spatial gradients, the Lagrangian averaged Hessian is 
a challenging object to compute. Several attempts have been made to model the Lagrangian 
averaged pressure Hessian in (5.2) by introducing a constitutive closure for it. This idea goes 
back to Leorat [24], Vieillefosse [25] and Cantwell [26] who assumed that the Eulerian pressure 
Hessian P is isotropic : see also [27, 28]. This assumption results in the restricted Euler equations 
(5.1) and (5.2) with 

P = -^tr(M 2 ), (5.3) 

where ld a b = 5 a b an d tr (Id) = 3 in three dimensions, so that taking the trace satisfies the relation 
required for incompressibility [30]. Conversely, one may assume that the Lagrangian pressure 
Hessian is isotropic. The latter assumption underlies the mean flow features of the tetrad model 
of Chertkov, Pumir and Shraiman [31]; see Chevillard and Meneveau [32] for a recent review of 
this approach and results on its use in turbulence diagnostics. 

A more general model that encompasses the mean flow features of both the restricted Euler 
equations and the tetrad model emerges from the transformation properties of the pressure 
Hessian under the Lagrangian flow map. The pressure Hessian transforms from Eulerian to 
Lagrangian coordinates by the flow map <pt : X — ► x(t) in which x(t) denotes the present 
position of a certain fluid particle, at time t, that started initially at position X = x(0) at time 
t = 0. This flow map preserves volume. Consequently, it is invertible and its Jacobian, the 
deformation gradient tensor, D l A (X,t) = dx l /dX A has unit determinant det(D) = 1. If it were 
frozen into the flow, the pressure Hessian would transform under the flow map as 



dx*dxJ dX%{t) ® dxJ(t) = ^ ° \ dX*dXB dX ® dX ) ■ (5 ' 4) 

This is the way that a Riemannian metric G transforms under a time-dependent change of spatial 
coordinates. Namely, the transformation of a Riemannian metric gives the evolving quantity 
G(t) in terms of the reference metric G(0) and the change of basis governed by the evolution of 
Jacobian matrix (D -1 )^(t) = dX A /dx l of the inverse flow map in tensor index notation as 

G ij (t) = G AB (0)(D- 1 )t(t)(D- 1 )f(t). (5.5) 

In short, a Riemannian metric in the Lagrangian reference configuration transforms under the 
flow map 4> t as 

G(t) = D- 1 (t)G(0)D- 1 (t), (5.6) 

where the material time derivative of D~ 1 (t)dx(t) = dx(0) determines the evolution of D -1 (t) 2 
For the pressure Hessian to transform as a Riemannian metric and also to satisfy the Poisson 
equation for its trace, it must take the following algebraic form 

P = -7^ tr ( M2 )- (5-7) 

Setting G(t) = Id recovers the restricted Euler equations of [25, 26], while setting G(0) = Id 
reformulates the mean flow part of the tetrad model [31]. 



2 The metric G(t) is called the Finger tensor in nonlinear elasticity. For its history and a modern application 
of the Finger tensor, see [33, 34, 29]. 
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In fact, the Poisson equation for pressure is satisfied for any choice of the nonsingular sym- 
metric matrix G = G T in equation (5.7) for the Hessian. Moreover, it may even be satisfied by 
choosing a linear combination of symmetric matrices in the form 



so long as an evolutionary flow law is provided for each of the symmetric tensors = Gj with 
(3 = 1, ... , N. Any choice of these flow laws would also determine the evolution of the driving 
term q& in the Riccati equation (3.13). 

6 Conclusions 

The review of rigid body rotations in §2 shows that, if handled properly, quaternions have 
a computational advantage over Euler angle formulations. For example, the Euler-Rodrigues 
formula (2.6) arises from a simple multiplication of quaternions. When this approach is applied 
to Lagrangian evolution equations, as in §3, it demonstrates that quaternions are a natural 
and efficient way of calculating the orientation and angular velocity of Lagrangian particles in 
motion through the introduction the concept of ortho-normal quaternion-frames travelling with 
each particle. For any problem in this class, knowledge of the quartet of 3-vectors (u, w, a, b) is 
sufficient for the application of Theorem 1, which is the paper's main result. The complexity of 
the various versions of the 3D Euler equations comes through the ortho-normal dynamics via the 
pressure field which is itself coupled back through Ap = —UijUj^. This has been discussed more 
fully in [35, 36]. Adaptations of these ideas when more physics is added to the Euler equations 
has been discussed in section §4. For the 3D Euler equations the pressure affects the coefficient 
Cb in the Darboux vector T> a b of Theorem 1 through the pressure Hessian P. The pressure 
Hessian also figures prominently in Lagrangian averaged models of velocity gradient dynamics 
in turbulence. These models employ approximate constitutive relations for the pressure Hessian. 
The potential applicability of the present quaternionic frame description in these Lagrangian 
averaged turbulence models was demonstrated in §5. 

It is also possible that the general formulation could be modified to include viscous effects, 
particularly if experimental data becomes available in a quaternionic format: the reader is 
referred to the review [37] and also to [38]. While the current formulation depends only on u> 
and not Vu>, the latter dependence could not be avoided if viscosity were included. An equivalent 
formulation for the compressible Euler equations ([39, 40]) may give a clue to the nature of the 
incompressible limit. 
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